library(sf) # working with spatial vector data
library(terra) # working with spatial raster data
library(tmap) # plotting maps
# load raster data
# the pipe operator |> is available for R version 4.1 or higher
rasters <- system.file("extdata/au/", package = "blockCV") |>
list.files(full.names = TRUE) |>
terra::rast()
install.packages("blockcv")
install.packages("blockCV")
# load raster data
# the pipe operator |> is available for R version 4.1 or higher
rasters <- system.file("extdata/au/", package = "blockCV") |>
list.files(full.names = TRUE) |>
terra::rast()
points <- read.csv(system.file("extdata/", "species.csv", package = "blockCV"))
head(points)
pa_data <- sf::st_as_sf(points, coords = c("x", "y"), crs = 7845)
tm_shape(rasters[[1]]) +
tm_raster(
col.scale = tm_scale_continuous(values = gray.colors(10)),
col.legend = tm_legend_hide()
) +
tm_shape(pa_data) +
tm_dots(
fill = "occ",
fill.scale = tm_scale_categorical(),
size = 0.5,
fill_alpha = 0.5
)
library(blockCV)
scv1 <- cv_spatial(
x = pa_data,
column = "occ", # the response column (binary or multi-class)
r = rasters,
k = 5, # number of folds
size = 360000, # size of the blocks in metres
selection = "random", # random blocks-to-fold
iteration = 50, # find evenly dispersed folds
progress = FALSE, # turn off progress bar
biomod2 = TRUE, # also create folds for biomod2
raster_colors = terrain.colors(10, rev = TRUE) # options from cv_plot for a better colour contrast
)
range <- cv_spatial_autocor(
x = pa_data, # species data
column = "occ", # column storing presence-absence records (0s and 1s)
plot = FALSE
)
range$range
scv2 <- cv_nndm(
x = pa_data,
column = "occ",
r = rasters,
size = 360000, # range of spatial autocorrelation
num_sample = 10000, # number of samples of prediction points
sampling = "regular", # sampling methods; it can be random as well
min_train = 0.1, # minimum portion to keep in each train fold
plot = TRUE
)
source("~/Dropbox/work/reimbursement and travel/202601 Qingdao/teaching/Topic 7/entropyweight.R")
dt <- read.csv("urbanvars.csv")
setwd(''/Users/268222h/Dropbox/work/reimbursement and travel/202601 Qingdao/teaching/Topic 7'')
setwd('/Users/268222h/Dropbox/work/reimbursement and travel/202601 Qingdao/teaching/Topic 7')
dt <- read.csv("urbanvars.csv")
View(dt)
dt <- read.csv("urbanvars.csv")
vars1 <- dt[, 2:8]
vars2 <- dt[, 9:10]
source("~/entropyweight.R")
source("~/Dropbox/work/reimbursement and travel/202601 Qingdao/teaching/Topic 7/entropyweight.R")
normlize.vars1 <- normalize(vars1)
normlize.vars1 <- normalize(vars1, direction = 1)
normlize.vars2 <- normalize(vars1, direction = 1)
normlize.vars2 <- normalize(vars2, direction = 1)
w1 <- entropyweight(normlize.vars1)
w1
w2 <- entropyweight(normlize.vars2)
w2
z1 <- colSums(t(normlize.vars1) * w1)
z2 <- colSums(t(normlize.vars2) * w2)
w1 <- entropyweight(normlize.vars1)
w1
w2 <- entropyweight(normlize.vars2)
w2
z1 <- colSums(t(normlize.vars1) * w1)
z2 <- colSums(t(normlize.vars2) * w2)
vars1 <- dt[, 2:8]
vars2 <- dt[, 9:10]
View(vars1)
normlize.vars1 <- normalize(vars1, direction = 1)
normlize.vars2 <- normalize(vars2, direction = 1)
w1 <- entropyweight(normlize.vars1)
w1
w2 <- entropyweight(normlize.vars2)
w2
z1 <- colSums(t(normlize.vars1) * w1)
z2 <- colSums(t(normlize.vars2) * w2)
View(normlize.vars1)
View(dt)
dt <- read.csv("urbanvars.csv")
vars1 <- dt[, 2:8]
vars2 <- dt[, 9:10]
normlize.vars1 <- normalize(vars1, direction = 1)
normlize.vars2 <- normalize(vars2, direction = 1)
w1 <- entropyweight(normlize.vars1)
w1
w2 <- entropyweight(normlize.vars2)
w2
View(normlize.vars1)
View(normlize.vars1)
View(normlize.vars2)
normlize.vars1 <- normalize(vars1, direction = rep(1,7))
normlize.vars2 <- normalize(vars2, direction = rep(1,2))
w1 <- entropyweight(normlize.vars1)
w1
w2 <- entropyweight(normlize.vars2)
w2
z1 <- colSums(t(normlize.vars1) * w1)
z2 <- colSums(t(normlize.vars2) * w2)
urbanindex <- (z1 + z2)/2
urbanindex <- data.frame("PID" = dt$PID, urban)
urban <- (z1 + z2)/2
urbanindex <- data.frame("PID" = dt$PID, urban)
View(urbanindex)
# Read shapefile and standardize PID type
sa <- st_read("studyarea.shp", quiet = TRUE) %>%
mutate(PID = as.character(PID))
library(sf)
library(dplyr)
library(ggplot2)
library(patchwork)
library(scales)
# Read shapefile and standardize PID type
sa <- st_read("studyarea.shp", quiet = TRUE) %>%
mutate(PID = as.character(PID))
setwd('/Users/268222h/Dropbox/work/reimbursement and travel/202601 Qingdao/teaching/Topic 7/generate urban index')
# Read shapefile and standardize PID type
sa <- st_read("studyarea.shp", quiet = TRUE) %>%
mutate(PID = as.character(PID))
# Standardize PID type in table
sp.urbanindex <- urbanindex %>%
mutate(PID = as.character(PID))
# Join attributes to polygons
sa2 <- sa %>%
left_join(sp.urbanindex %>% select(PID, urban), by = "PID")
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_viridis_c(option = "turbo") +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_viridis_c(option = "magma") +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_viridis_c(option = "viridis") +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_viridis_c(option = "rocket") +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_viridis_c(option = "mako") +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_viridis_c(option = "inferno") +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_gradientn(
colours = c("#2c7bb6", "#00a65a", "#d7191c")
) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_gradient2(
low = "blue",
mid = "green",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE)
) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_gradient2(
low = "blue",
mid = "orange",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE)
) +
theme_bw()
gplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_gradient2(
low = "blue",
mid = "orange",
high = "red",
lowpoint = min(sa2$urban),
midpoint = median(sa2$urban),
highpoint = max(sa2$urban)
) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray") +
scale_fill_gradient2(
low = "blue",
mid = "orange",
high = "red",
lowpoint = min(sa2$urban),
midpoint = median(sa2$urban),
highpoint = max(sa2$urban)
) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "blue",
mid = "orange",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
oob = squish,
na.value = "grey90",
name = "Urban index"
) +
coord_sf(datum = NA) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "blue",
mid = "white",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
oob = squish,
na.value = "grey90",
name = "Urban index"
) +
coord_sf(datum = NA) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "blue",
mid = "white",
high = "darkred",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
oob = squish,
na.value = "grey90",
name = "Urban index"
) +
coord_sf(datum = NA) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "darkblue",
mid = "white",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
oob = squish,
na.value = "grey90",
name = "Urban index"
) +
coord_sf(datum = NA) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "blue",
mid = "white",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
oob = squish,
na.value = "grey90",
name = "Urban index"
) +
coord_sf(datum = NA) +
theme_bw()
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "blue",
mid = "white",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
na.value = "grey90"
) +
theme_bw()
###############################################################################
# 9) Visualise spatial distribution of the urbanisation index
#    A diverging colour scale (blue–white–red) highlights low, medium,
#    and high levels of urbanisation
###############################################################################
ggplot(sa2) +
geom_sf(aes(fill = urban), color = "gray60", linewidth = 0.1) +
scale_fill_gradient2(
low = "blue",
mid = "white",
high = "red",
midpoint = median(sa2$urban, na.rm = TRUE),
limits = range(sa2$urban, na.rm = TRUE),
na.value = "grey90"
) +
theme_bw()
